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ABSTRACT 

Certain algebraic combinations of single scattering albedo and solar radiation reflected from, or transmitted 
through, vegetation canopies do not vary with wavelength. These “spectrally invariant relationships” are the 
consequence of wavelength independence of the extinction coefficient and scattering phase function in veg- 
etation. In general, this wavelength independence does not hold in the atmosphere, but in cloud-dominated 
atmospheres the total extinction and total scattering phase function vary only weakly with wavelength. This 
paper identifies the atmospheric conditions under which the spectrally invariant approximation can accu- 
rately describe the extinction and scattering properties of cloudy atmospheres. The validity of the as- 
sumptions and the accuracy of the approximation are tested with ID radiative transfer calculations using 
publicly available radiative transfer models: Discrete Ordinate Radiative Transfer (DISORT) and Santa 
Barbara DISORT Atmospheric Radiative Transfer (SBDART). It is shown for cloudy atmospheres with 
cloud optical depth above 3, and for spectral intervals that exclude strong water vapor absorption, that the 
spectrally invariant relationships found in vegetation canopy radiative transfer are valid to better than 5%. 
The physics behind this phenomenon, its mathematical basis, and possible applications to remote sensing and 
climate are discussed. 


1. Introduction 

Recently several papers reported the discovery of 
spectrally invariant behavior in some simple algebraic 
combinations, called “spectral invariants,” of single 
scattering albedo and solar radiation reflected from or 
transmitted through vegetation canopies (Knyazikhin 
et al. 1998, 2005; Huang et al. 2007). The spectral in- 
variant phenomenon is clearly seen in radiative mea- 
surements and remote sensing data (Panferov et al. 2001; 
Wang et al. 2003). The phenomenon was theoretically 
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explained and numerically simulated using radiative 
transfer theory (Knyazikhin et al. 2011; Smolander and 
Stenberg 2005). There are three key parameters that 
characterize the radiative transfer process: the extinction 
coefficient, scattering phase function, and single scatter- 
ing albedo. In vegetation, the optical distance between 
two points within the canopy does not depend on wave- 
length because the scattering elements are much larger 
than the wavelength of solar radiation (Ross 1981). And 
the canopy scattering phase function is also wavelength 
independent since it is determined by large scattering 
elements such as twigs and leaves. Thus, of the three key 
variables, the single scattering albedo is the only one 
with significant wavelength dependency. This allows for 
a natural separation between structural and spectral pa- 
rameters of radiative transfer: the extinction coefficient 
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and scattering phase function are purely structural, while 
the single scattering albedo is purely spectral. 

This separation does not hold for the atmosphere 
where extinction and scattering due to Rayleigh mole- 
cules and aerosols, as well as gaseous absorption, depend 
strongly on wavelength. However, under some cloudy 
atmospheric conditions the physical processes that vary 
weakly with wavelength (e.g., cloud scattering) may dom- 
inate. Then, for a large range of wavelengths, the extinc- 
tion coefficient and the scattering phase function may 
be weakly spectrally variable. 

In loose analogy with the structural/spectral separa- 
tion that occurs in vegetation radiative transfer, the 
standard method for cloud remote sensing from space 
(Nakajima and King 1990) retrieves optical depth and 
droplet effective radius. Optical depth is the structural 
parameter with no significant wavelength dependence, 
while effective radius strongly affects single scattering 
albedo (Twomey and Bohren 1980). The Nakajima-King 
method decomposes the radiative measurements into 
structural (optical depth) and spectral (effective radius) 
parameters. It takes advantage of the fact that radiation 
in the visible, where particles do not absorb, is much more 
sensitive to cloud optical depth than to droplet size. In 
contrast, in the near infrared, where water absorbs, the 
radiance changes substantially with effective radius and 
is insensitive to cloud optical depth, at least for thicker 
clouds. 

The above analogy between cloud and vegetation 
remote sensing is not perfect because generally cloud 
droplet size also affects the scattering phase function. 
However, we can still ask how well the spectral invariant 
approach in vegetation canopies works for cloudy atmo- 
spheres. The goal of this paper is to answer that question. 

The plan of the paper is as follows. In the next section 
we sketch the general concept of spectral invariance (SI) 
and demonstrate its validity with simple Discrete Ordi- 
nate Radiative Transfer model (DISORT) calculations. 
Section 3 looks at the extent to which the spectral vari- 
ability of the extinction and scattering properties in cloudy 
atmospheres meets the needs of spectral invariance the- 
ory. Section 4 illustrates the spectrally invariant behavior 
of radiances and fluxes and discusses its physical in- 
terpretation. Accuracy of the spectral invariant approxi- 
mation is discussed in section 5, and section 6 briefly 
describes its possible implementations in remote sensing 
and climate modeling. Finally, the appendix provides 
a physical and mathematical basis for spectral invariance. 

2. Spectral invariance 

Spectral invariance states that simple algebraic com- 
binations of the spectra of single scattering albedo and 


f ^ ^ incident beam 






r \ p ' O 

recollision W 


%P&) 





recollision 

/ 

T J 


non-reflccting surface 


FIG. 1. Schematic of radiative transfer process. The term T dir is 
the fraction of photons in the incident beam that reach the surface 
without interacting. A fraction fo = 1 - T dit therefore interacts with 
the medium. With probability w 0 a. these photons are scattered and 
then either interact again (with probability p) or escape the me- 
dium in direction ft [with probability p(ft)]. 


reflected and/or transmitted radiation eliminate the de- 
pendences on wavelength (Huang et al. 2007; Knyazikhin 
et al. 2011). In many cases spectral dependency can be 
compressed into a linear relationship 

y(A) = ax(\) + b, (1) 

where parameters a and b are independent of wave- 
length A, and functions y(A) and x(\) are algebraic com- 
binations of measured spectra. 

Spectrally invariant relationships were originally de- 
rived for vegetation canopies with “black” soil. Here we 
briefly summarize the “nontraditional” derivation (e.g., 
Schull et al. 2007) of radiative transfer in an absorbing 
and reflecting layer, following the diagram in Fig. 1. 
We should emphasize that we assume that the standard 
diffuse-direct transformation of the radiative transfer 
equation has been applied; we are only dealing with the 
diffuse part. 

Let Tdir be the direct transmittance — that is, the frac- 
tion of photons from the incident beam that pass through 
the layer in Fig. 1 without colliding. With probability 
cuq a collided photon will be scattered and will either 
recollide or escape the layer in a given direction fl with 
probabilities p and p(fl), respectively. Here oj 0 , p, and 
p(H) are the single scattering albedo, recollision proba- 
bility, and directional escape probability, respectively. 

While oi 0 is a well-known parameter in atmospheric 
radiation, p and p(Cl) are less known and thus require 
some explanation: p is the conditional probability that 
a scattered photon will interact with the medium again 
(Smolander and Stenberg 2005), and p(il) is the condi- 
tional probability that a scattered photon will leave the 
medium in the direction H. The parameters p and p are 
simply related. Let us assume for simplicity that T dir = 0. 
If p is the probability of recollision, then 1 - p is the total 
escape probability and (e.g., Schull et al. 2011) 
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where p, is the cosine of a polar angle of Cl. Here p(f 1)/ 
(1 - p) is the probability density of escape in fl as a re- 
sult of one interaction; it can be considered as the whole 
domain counterpart of the scattering phase function. In 
general, p and p depend on the successive order of scat- 
tering; however, this dependence weakens as the number 
of interactions increases (Huang et al. 2007). 

If we assume that p and p do not depend on the 
scattering order, the probability that a photon will es- 
cape the layer after m interactions is i Q pp m ~ l <o (ft, where 
io = 1 - T t jjr is the fraction of incoming photons that 
collide within the layer (Fig. 1). For a purely absorbing 
surface (zero lower boundary condition), the fraction of 
exiting photons in a given Cl is simply the sum of these 
probabilities over scattering order m; that is. 


= P(°) W 0A*0 + P(tt)6&p/ 0 + 
+ + .. 


0A 


1 - pw 0A 


p(Cl)i Q . 


( 3 ) 


So far we have suppressed wavelength dependency in all 
quantities except / and w 0 ; w 0 is assumed to be the only 
spectrally varying parameter. Equation (3) can be re- 
arranged as 


^ = p/ A (ft) + p(ft), 0 . (4) 

"0A 

If p and pi 0 do not depend on wavelength, Eq. (4) 
qualifies as an SI relationship [Eq. (1)]. 

Under what conditions is Eq. (4) valid in the atmo- 
sphere? In vegetation canopies the extinction coefficient 
a is wavelength independent — that is, 

o-(A) = a (5a) 

— and the scattering phase function P does not depend 
on wavelength either: 




(b) I(x,SZA,VZA,VAA) 

Fig. 2. Linear relationship between the ratio Ik r) 0 and / (/ is re- 
flected radiance) for r = 10, and for seven values of <u 0 from 0.8 to 
1.0. A Henyey-Greenstein scattering phase function with asym- 
metry factor g = 0.85 is used. The surface is black. The ID radiative 
transfer code DISORT is used to calculate values of /(r, SZA, 
VZA, VAA) where SZA, VZA, and VAA are the solar zenith, 
viewing zenith, and viewing azimuth angles, respectively. SZA is 
60°. (a) VZA = 0°. (b) VZA = 60°; VAA = 0° (black circles) and 
VAA = 180“ (gray squares). 


E(n,n\A) = P(Cl,Cl'). (5b) 

The single scattering albedo wq* thus becomes the only 
spectrally varying parameter in the radiative transfer equa- 
tion. As a result, p and p(Cl) become wavelength in- 
dependent (Knyazikhin et al. 2011). In this case, Eq. (4) 
forms a straight line on the / A /eu 0A versus / A plane where 
slope and intercept give p and escape factor pi 0 . (A for- 
mal definition of the recollision probability in terms of the 
3D radiative transfer equation is given in the appendix). 


Figure 2 illustrates the validity of Eq. (4) using the ID 
radiative transfer code DISORT (Stamnes et al. 1988) 
with the Henyey-Greenstein scattering phase function 
for a cloud with optical depth r = 10 and different values 
of w 0 between 0.8 and 1. Although the slope a (an ap- 
proximation of p) is mostly determined by the medium 
internal structure given by a and P (Smolander and 
Stenberg 2005), it also depends weakly on the solar- 
viewing geometry as shown in Fig. 2b. However, the 
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viewing geometry mostly affects the intercept of the lin- 
ear relationship b, an approximation of p(fl). 

Obviously, the assumptions in Eqs. (5a) and (5b) are 
not met for the atmosphere where both the extinction 
due to molecular and aerosol scattering and the gaseous 
absorption strongly depend on wavelength. However, 
for cloudy conditions and certain wavelength ranges, 
the physical processes that vary only weakly with A may 
dominate. As a result, the atmospheric extinction co- 
efficient and the scattering phase function may be ap- 
proximately constant there. Section 3 explores this 
possibility. 

3. Spectral variability of the extinction and 
scattering properties in real atmosphere 

Let us calculate the total extinction coefficient as a 
function of wavelength: 

N N 

°a = = + 

k ~ 1 k = 1 

(6a) 

Here index k refers to the fcth atmospheric constituent, 
s to scattering, and a to absorption; for example, Rayleigh 
molecules (cr SjA ^ 0, tj aX = 0), aerosol particles (cr SiA / 0, 
o aA ^ 0), cloud droplets (<t vA ^ 0, (j aA ^ 0), and gaseous 
absorbers (a sA = 0, <r a>A ^ 0). In the radiative transfer 
equation the total (or effective) single scattering albedo 
is defined as 

1 n ^ N ^ 

"OA = = = ( 6b ) 

°A k=l °A k= 1 (r x 

For simplicity, we will characterize the scattering phase 
function with the asymmetry parameter: 

N 

( 6c > 

Instead of the total extinction coefficient we will use the 
total optical depth: 

N S 

t a = 2 r u = X (t jU + r flU ) - + r flA . (6d) 

/c=l k- 1 

For three values of cloud optical depth, Fig. 3 illustrates 
the total optical depth r A = trx + t ax + r CA + t ga and 
each of its components: Rayleigh molecules trx, aero- 
sols tax, clouds t ca , and gas t ca . (Since r CA only weakly 


varies with A , we will drop the subscript A , understanding 
that tc = tc, a =o, 55 M m)- We see that for a thin cloud 
(Fig. 3b) spectrally variable aerosol, air molecules, and 
gas dominate in the total optical depth and there are few 
spectral intervals where t a does not vary over a large 
range. However, for thicker clouds (r c ^ 3), a spectrally 
“flat” cloud dominates, making the total optical depth 
almost insensitive to spectral variations in r RX and t ax 
outside the water vapor absorbing bands. If strongly water 
vapor absorbing wavelengths are excluded, the extinction 
coefficient becomes weakly dependent on wavelength 
and the first condition for the applicability of SI stated 
in Eq. (5a) is approximately met. 

The total asymmetTy parameter of the scattering phase 
function (Fig. 4a) for optically thin clouds is strongly 
affected by molecular scattering, especially for shorter 
wavelengths. For thicker clouds, it exhibits much weaker 
variability with wavelength, suggesting that the second 
condition for SI, Eq. (5b), is also approximately met. 

If the total extinction coefficient and the scattering 
phase function are only weakly dependent on wave- 
length, the spectral scattering and absorption processes 
in clouds are almost entirely determined by spectral 
variation in o)qx defined by Eq. (6b). Figure 4b shows 
coox for three cloud optical depths: tc = 0.5, 3, and 10. 
Note that, in contrast to the cloud liquid water single 
scattering albedo that varies spectrally between 1 and 
0.96 depending on droplet size, the total single scattering 
albedo is dominated by gaseous absorbing spectral bands 
and may fluctuate between 0 (gas only) and 1 (no ab- 
sorption), especially for thinner clouds. This also sug- 
gests that while thick clouds suppress spectral variability 
in the total optical depth and scattering phase function, 
w 0A varies significantly because of its strong sensitivity to 
the absorbing and scattering properties of the atmo- 
spheric components. This feature makes it possible to 
implement the SI approach by directly relating radiative 
spectral properties of the atmosphere to the spectrum 

of £t>0A- 

To increase the range of w 0A while keeping low spec- 
tral variation in r A and g A , we will begin our analysis with 
only those wavelengths that have wqa > 0.8. Table 1 
summarizes spectral variation in t a and g A for those 
wavelengths. We can see that for t c s 3, clouds suppress 
spectral sensitivity of r A to trx and t ax optical depths. As 
a result, for wavelengths with co 0A s: 0.8, the spectral 
variation in r A and g A does not exceed 5% and 2%, re- 
spectively. 

Finally we recall that Eq. (3) was derived under the 
assumption that the recollision and directional escape 
probabilities do not depend on the scattering order (see 
Fig. 1). In a general case, the spectral invariants vary 
with the number of successive collisions; however, they 
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Fig. 3. Spectrum of total optical depth of a cloudy atmosphere as defined by Eq. (6d), ac- 
counting for air molecules (Rayleigh scattering), cloud droplets (scattering and absorption), 
aerosol particles (scattering and absorption), and gases (absorption). Optical depth of aerosol at 
A = 0.55 jum is t a = 02. (a) Three values of cloud optical depth: r c = 0.5, 3, and 10 at A =■ 
0.55 /xm. (b) Zoom of spectra for r c = 0.5, cutting off at r = 2 as indicated by the dotted line in (a). 


asymptote quickly as the number of collisions increases 
(Huang et al. 2007). For conservative scattering, Davis 
and Marshak (1997) estimated the number of forward 
scatterings required for particle to lose all memory of its 
original direction of travel as (1 - g)~ ] ~ 6-7. 

4. Radiative transfer calculations 

To compute quantity / A (0) in Eq- (3) and its hemi- 
spherically integrated counterparts, the reflectance R x 
and transmittance T x , for wavelength range 0.4-2.5 fim 
with 10-nm resolution, we use the publicly available 


Santa Barbara DISORT Atmospheric Radiative Trans- 
fer (SBDART) code (Ricchiazzi et al. 1998) with a solar 
zenith angle of 45°. The other key parameters used in our 
simulations are listed in Table 2 and their choice is dis- 
cussed in Chiu et al. (2010). 

Let us come back to Eq. (3), multiply both sides by the 
cosine of a polar angle, and integrate over all 47 t solid 
angle. Accounting for Eq. (2), we get 

w _ 0 - ~ pKa 

A 1 “ P ™ oa ’ 


( 7 ) 
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0.5 1 1.5 2 2.5 

wavelength, X (pm) 

(b) 


FIG. 4. (a) Spectra of g A for three cloud optical depths: r c = 0.5, 3, and 10 with T/ti = n = 
0.2. (b) Spectra of wqa for the same three cases as in (a). See Eqs. (6b) and (6c) for definitions of 
'‘total” in each case. 


where 

W A d T dit ) =R, + (T k - T dir ) 


(8a) 


corresponds to the single scattering co albedo ^oa = 1 — 
iu 0 A- From Eq. (7) we can prove that 


is the total radiation that has been reflected R A from, or 
diffusively transmitted (7\ - T diT ) through, a cloudy 
layer. For sufficiently thick clouds the direct radiation 
T dir « 0 and Eq. (8a) becomes 

W'a = + A- ( 8fe ) 

Here W A is the whole domain counterpart of wqa while 


because p < 1, and thus 

Aa “ fl 0A' 

Interestingly, the /l 0 A-t°- a 0 A ratio estimates the average 
number of photon interactions as a function of wave- 
length (Panferov et al. 2001; Knyazikhin et al. 2005): 


Aa 

_ 1 

"A 

a 0A 

“ 1 

3 

3 

1 



Aa = i - A 


* 1. (9) 
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TABLE 1. Variations in the total optical depth and scattering phase function asymmetry parameter corresponding to wavelengths for 
which the single scattering albedo exceeds 0.8. The second row is for wavelengths between 0.7 and 2.5 /xm. In the “x ± y" notation in 
columns 3 and 5, x and y stand for the mean and the standard deviation, respectively. 


Cloud optical 
depth 

No. of wavelengths 
from 0.4 to 2.5 /xm 
(cut of total 210) 

Total optical 
depth 

Std dev to mean ratio for 
total optical depth (%) 

Total asymmetry 
parameter 

Std dev to mean ratio 
for total asymmetry 
parameter (%) 

0.5 

97 

0.72 ± 0.12 

17 

0.79 ± 0.07 

9 

0.5 

67 

0.65 ± 0.05 

8 

. 0.83 ± 0.02 

2 

3 

154 

3.34 ± 0.16 

5 

0.85 ± 0.02 

2 

5 

166 

5.44 ± 0.24 

4 

0.86 ± 0.01 

1 

10 

179 

10.66 ± 0.44 

4 

0.86 ± 0.01 

1 


The final equality in Eq. (9) follows directly from Eq. (7). 

Note that the reciprocal of n A is linearly related to a> 0A . In 
addition, the product of and the photon mean free path 
l/cr estimates the total photon path length: or 


L 


a 


”A _ 1 * 

O- 0-1 - P<U 0A ' 


( 10 ) 


—L^aW^+b (11a) 

w o\ 


W = 

A 1 


bo), 


0A 


aco, 


0A 


(Hb) 


Equation (10) shows the explicit dependence on p and 
allows us to estimate L A using explicit numerical models 
like SBDART [or the spherical harmonics discrete or- 
dinate method (SHDOM); Evans 1998] rather than sta- 
tistical models like Monte Carlo. 

Figure 5 shows spectral variations of n A and W A . There 
are two remarkable points here. First, in the spectral 
interval 0.4-0.55 pm, n A increases from 17 to 23 for 
a thick cloud ( tc =10) but slowly decreases for a thin 
one (rc = 0.5). This behavior follows from the fact that 
for thicker clouds w 0 a (Fig. 4b) decreases much more 
slowly than W A [see Fig. 5b and Eq. (10)], while for thin 
clouds, because of absorbing aerosols, w 0 a decreases 
faster than W\. The recollision probability controls the 
response rate of VV A to variation in woa- "The second re- 
markable point about Fig. 5 is the negligible effect 
of drop size except at wavelengths where liquid water 
absorbs strongly. In the spectral interval 2. 1-2.2 /urn, 
\6-pm droplets have an average of 13 photon inter- 
actions for tc — 10, while for 4-/xm (and therefore less 
absorbing) droplets, n A goes up to 18. Finally, for the 
strongly water vapor-absorbing spectral bands around 
1.4 and 1.8 pm, W A = 0 and n A approaches 1. 

Figure 6 plots points only for wavelengths for which 
cdoa > 0.8 (see Table 1); there are 67 (r c = 0.5), 154 (t c = 
3.0), and 179 (t c = 10.0) such wavelength intervals. 
Figure 6a is similar to Fig. 5b. Figure 6b shows the ratio 
WJat ox plotted against W x , where the slope increases as 
we go from aerosol- and molecular-dominated condi- 
tions ( tc = 0.5) toward cloud-dominated conditions 
(tc = 3, 10). The linear fit is almost perfect for t c = 10. 
This confirms that SI holds for the cloud dominated 
cases — that is, 


Note that the sum of slopes a and intercepts b in Fig. 6b 
is almost exactly equal to unity. In the zero-absorption 
limit where the total u/oa = 1, W A = 1 and consequently 
a -f b — 1 from Eq. (11a). 

Note that Eq. (lib) coincides with Eq. (7) if a = p and 
b = 1— p. Thus Fig. 6b numerically confirms the validity 
of Eq. (7). Finally, Eq. (11a) is equivalent to 


Table 2. SBDART parameters. 


Parameters 

Values used in model 

Spectral 

Lower wavelength limit 

0.4 pm 

Upper wavelength limit 

2.5 pm 

Spectral resolution 

0.01 /xm 

Solar 

Solar spectrum 

MODTRAN_3 

Solar zenith angle 

45° 

Atmosphere 

Atmospheric profile 

Midlatitude summer 

Integrated water vapor amount 

3 cm 

Integrated ozone concentration 

0.324 atm-cm 

Surface 

Surface albedo 

0 

Cloud 

Cloud-base altitude 

1 km 

Cloud optical depth at 0.55 pm 

0.0, 0.5, 3, 5, 10 

Cloud drop effective radius 

4 and 16 pm 

Cloud phase function 

Henyey-Greenstein 

Aerosol 

Aerosol type 

Rural 

Aerosol optical depth at 0.55 /xm 

0.2 

Aerosol phase function 

Henyey-Greenstein 

Visibility at 0.55 pm 

23 km 

Relative humidity used in the 

80% 

boundary layer aerosol model 
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0.5 1 1.5 2 2.5 

- . wavelength, X (pm) 
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Fig. 5. (a) Spectra of estimates of n A for clouds with r c = 0.5, 3, and 10 (r eff = 4 and 16 /*m); 
Ta . a = 0 . 55^01 = 0-2- Dashed horizontal line corresponds to one interaction, (b) Spectra of W A for 
the same cases as in (a). 


W 

-±= p w k + (l- p ). ( 12 ) 

"(U 

Based on Eq. (12) the following interpretation of the 
linear relationship observed in Fig. 6b can be given. The 
terms wqa and W A quantify the scattering events result- 
ing from one and multiple interactions, respectively. 
Their ratio is the probability that a scattered photon will 
escape the layer as a result of multiple interactions. Part 
of the single scattered photons will escape without fur- 
ther interactions with the probability (1 - p), and the 
rest with the probability pW A . 

Interestingly, while both W A and o>oa depend on 
droplet radius (at least for the 2.1-2. 2- pm spectral 


region, as seen in Fig. 5b), the slopes a (or p) are almost 
independent of droplet radius. Hence, the increase in n A 
for smaller droplets is due entirely to larger £u ()A , as fol- 
lows from Eq. (9). 

Finally, Fig. 7 shows the same kind of plot as in Fig. 
6b, but for zenith and nadir radiances instead of W A . 
We see that (i) again, the behavior is very close to 
linear, at least for the cloud-dominated cases; (ii) for 
r c = 10 the slopes for W A and for the two radiances are 
almost identical, while for the other cases the slopes are 
different; and (iii) the sum of a and b is not necessarily 
equal to 1. 

In summary, for the cloud-dominated cases, both ra- 
diance and fluxes show evidence of SI as predicted by 
Eqs. (4) and (12). 
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FIG. 6. (a) Spectra of W A for clouds with Tc = 0.5, 3, and 10 (r eff = 16 /xm); T, 4iA =o.55 (i ra = 0.2. 
Doglegs occur because only wavelengths for which w 0A > 0.8 (see Table 1) are plotted, 
(b) Ratio W\l a»<u is plotted against W h for the wavelengths from (a). To reduce variability in r A 
and g A (see Table 1), wavelengths less than 0.7 /xm were not used for r c = 0.5. 


5. Accuracy of the spectrally invariant relationships 


To address the question of accuracy in the SI rela- 
tionships, we first introduce a more accurate wavelength- 
selection method than the criterion cu^ > 0.8 (Table 1) 
used in section 3. This method will be applied to obtain 
the slope and intercept. Second, we reconstruct the ra- 
diance spectra using the single scattering albedo and SI 
parameters and compare them with those calculated 
with the SBDART code. 

Solving Eq. (7) for p yields 


Px = 


OJ 0A 



(13) 


Our goal is to find the spectral intervals for which vari- 
ation in p A is small. Figure 8a shows variation in p A with 


wavelength. The sharp peaks around the absorbing 
bands are correlated with the total optical depth (cf. 
Fig. 3). In the spectral interval between 0.4 and 0.55 /xm, 
p A increases for thick clouds (re & 3) and decreases for 
a thin one (t c = 0.5). For A > 0.55 pm, p A becomes 
relatively flat outside the absorbing bands. Figure 8b 
shows the occurrences of the p A values from Fig. 8a while 
Tables 3 and 4 summarize their statistics for wave- 
lengths corresponding to the maxima in Fig. 8b. For 
such wavelengths, the W x /<o ox versus W x relationships 
are more linear than the / A /cu 0A versus / A relationships 
and their sums of slope and intercept are closer to unity. 
The tables show that the statistics of frequency peaks is 
only weakly sensitive to cloud drop size (as can be seen 
also in Fig. 8b). 

Equations (3) and (7) can be combined as 
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FIG. 7. As in Fig. fib, but for (a) / A (H) in zenith (ft = i ) direction 
and (b) / A (X1) in nadir (fl = T) direction, where l is reflected ra- 
diance. 


= K(n)Wi , (14a) 

L ''dir 

which, using Eq. (8a), is equivalent to 

/ A (o) = K (n)(R K + r A - r djr ), (14b) 

where *(fl) = p(fl)/(l - p) maps the total diffuse out- 
going flux into radiance [see Eq. (2)]. If the optical depth 
and the scattering phase function are independent of 
wavelength, the coefficient k does not depend on wave- 
length either and varies with fl (and the direction of 


incident radiation £L 0 ). Its physical interpretation is sim- 
ilar to the scattering phase function but formulated for 
the entire layer (i.e., the fraction k of scattered photons 
W x leaves the layer in a given fl). Figure 9 shows occur- 
rence of k values for nadir radiance. Unlike p values, k 
maxima become very sensitive to the drop sizes. If one 
selects only wavelengths that correspond to the maxima 
of both p and k, the differences in slopes between the 
two cases (hemispherical and directional) will be much 
smaller, as the theory predicts. 

Figure 10a illustrates the improvement of the line- 
arity in the WJo ) 0 A versus W A relationships with rejec- 
tion of wavelengths, where t a substantially fluctuates 
(see Fig. 3a). In addition to our standard case of w 0 a > 
0.8 with 179 wavelengths (Table 1), we also plotted the 
wqa thresholds of 0.9 with 156 wavelengths. As expected, 
the increased woa threshold slightly improves the re- 
gression coefficient and decreases the slope (thus in- 
creasing the intercept). 

The improved linearity makes the reconstruction of 
W x more accurate, especially for the shorter wavelengths 
(Fig. 10b). However, further increase in the threshold 
does not necessarily improve the accuracy. This is due to 
the reduced information on cloud liquid water absorp- 
tion by removing longer wavelengths (around 1.6 and, 
especially, 2.1 /im). Figure 11 illustrates the RMS and 
the bias 

RMS = M<W» - W IXaw f (15a) 

m m 

bias = — X W, A - -X W , . (15b) 

m-ti a m \ ti u - a PP r v ' 

for 10 approximations corresponding to the increased 
o)qx thresholds (m = 211). We see that there are few 
changes for the thresholds above o> 0 a = 0.97. This is be- 
cause most of the cloud liquid water absorbing wave- 
lengths around 2.1 pm have been removed and the range 
of juq* has been substantially decreased. The best ap- 
proximation is reached for the thresholds of 0.86 and 
0.92-93 for the bias and the RMS, respectively. Note 
that the threshold that gives the minimum RMS is very 
close to the one that is based on the maximum p values 
occurrence (see Table 4). 

What is the optimal number of wavelengths needed to 
obtain the recollision probability and what are those 
wavelengths? To increase the dynamic range of <o Q and 
thus the corresponding range of the I/coq versus I scat- 
terplot, one needs to use wavelengths with the largest 
range of cloud water absorption. In Fig. 12a we used 
three wavelengths: 0.4 pm (total o> 0 = 0.999, t = 10.57, 
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Fig. 8. (a) Spectra of p A for r c = 0.5, 3, and 10; r t « - 16 /im. (b) Histograms of p A from (a) for 
r ett = 4 and 16 with bin width 0.02. 



g = 0.835), 1.6 pun (w 0 = 0.982, r = 10.46, g = 0.861), and 
2.1 fj . m (wq = 0.948, r = 10.66, g = 0.866). The recon- 
structed W K for all wavelengths from 0.4 to 2.5 pan is 
shown in Fig. 12b. The RMS = 0.046 (~9% of the mean 


value of W x = 0.53) and the bias = —0.006. Note that 
accuracy of the approximation is not very sensitive to 
a selection of a particular wavelength in the neigh- 
borhoods of 0.4, 1.6, and 2.1 /am. Finally, it seems 


Table 3. Statistics of the ‘‘p A peak” values from Fig. 8b for r eU - = 4 p.m. Bin size is 0.02. The cloud-free column (t c = 0.0) is added for 
comparison. “Var” indicates coefficient of variation (the ratio of the standard deviation to the mean). 


Cloud optical depth (droplet size r eff = 4 pm) 



Parameters 

0.0 

0.5 

3.0 

5.0 

10.0 

Bin 


0.02-0.04 

0.36-0.38 

0.82-0.84 

0.90-0.92 

0.94-0.96 

No. of points 


27 

37 

59 

71 

115 

Mean p 


0.030 

0.371 

0.832 

0.910 

0.950 

Varp (%) 


17.0 

1.4 

0.7 

0.6 

0.6 

Min cjjw 


0.163 

0.722 

0.873 

0.884 

0.937 

Max woa 


0.655 

0.963 

0.998 

0.998 

0.999 

Var o>oa ( % ) 


47.4 

7.8 

3.9 

3.1 

1.4 

Min t a 


0.072 

0.630 

3.22 

5.27 

10.32 

Max r A 


0.170 

0.829 

3.90 

6.90 

12.51 

Var t a (%) 


24.7 

8.7 

5.8 

8.6 

5.5 

Correlation coefficient R 2 of linear relationship 

0.673 22 

0.998 63 

0.999 85 

0.999 89 

0.999 91 
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TABLE 4. As in Table 3, but for r eff = 16 j*m. “Var" indicates coefficient of variation (the ratio of the standard deviation to the mean). 


Parameters 


Cloud optical depth (droplet size r eff — 

16 Jim) 


0.0 

0.5 

3.0 

5.0 

10.0 

Bin 

0.02-0.04 

0.36-0.38 

0.82-0.84 

0.90-0.92 

0.94-0.96 

No. cf points 

27 

38 

65 

72 

104 

Mean p 

0.030 

0.367 

0.832 

0.909 

0.950 

Var p (%) 

17.0 

1.4 

0.7 

0.5 

0.6 

Min u 0A 

0.163 

0.683 

0.868 

0.869 

0.937 

Max wq* 

0.655 

0.964 

0.997 

0.998 

0.999 

Var ooa (%) 

47.4 

10.9 

3.7 

3.6 

1.5 

Min r A 

0.072 

0.587 

3.15 

5.21 

10.2 

Max t a 

0.170 

0.818 

3.16 

5.95 

11.0 

Var r A (%) 

24.7 

10.3 

3.8 

2.9 

1.5 

Rr of linear relationship 

0.673 22 

0.998 94 

0.999 86 

0.999 91 

0.999 91 


counterintuitive that the case with only three selected 
wavelengths (Fig. 12b) gives a better approximation than 
the one that uses more than 100 wavelengths (Fig. 11). 
However, the case with a large number of wavelengths is 
biased toward the visible spectral interval at the expense 
of near-infrared wavelengths. In addition, using a larger 
number of wavelengths leads to more violations of the 
assumptions stated in Eqs. (5a) and (5b). Thus more 
wavelengths do not necessarily provide a better approx- 
imation. 

6. Possible applications 

Guided by applications already being pursued in veg- 
etation canopy radiative transfer (Ganguly et al. 2008; 
Schull et al. 2011) we think the top applications are 

(i) testing the consistency of remote sensing retrievals, 

(ii) broadband calculations for climate models, (iii) fill- 
ing missing spectral data, and (iv) testing 3D radiative 
transfer codes. Below we briefly discuss all four. 


a. Testing the consistency of remote sensing retrievals 

Satellite remote sensing often provides particle ef- 
fective radius r eEf at three liquid water absorbing wave- 
lengths: 1.6, 2.1, and 3.7 /rm. From r e ff one can calculate 
the “effective” aj 1M for each wavelength. Now Eq. (4) 
can be used to validate the physical consistency of the 
retrieved r e tf values. Indeed, according to the SI theory, 
the ratio of / A to cooa should be approximately linear in 
Z A . Lack of linearity may call into question the correct- 
ness of retrievals. 

b. Broadband calculations for climate models 

For climate models it is assumed that if r ef f is known, 
the total w 0A can be calculated (or parameterized; see, 
e.g., Slingo 1989) as a function of wavelength (or spectral 
band). Using the assumptions of the SI approach, one 
can calculate fluxes for only a few (at least, two or three) 
wavelengths and obtain the two wavelength-independent 
variables (a and b). Fluxes for all other wavelengths could 



0.1 


0.2 


0.3 


0.4 


0.5 


0.6 


0.7 


K-value 

Fig. 9. As in Fig. 8b, but for k values. The term x(fl) = p(il)/(l — p) maps the total diffuse 
outgoing flux into radiance. Also, p and ?(Cl) are the recollision probability and the directional 
escape probability, respectively [see Eq. (2)]. 
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(b) wavelength, \(pm) 

FIG. 10. (a) Relationships of Wa/wq* vs W a for wavelengths such that > 0.8 (179 values) 
and taoA > 0.9 (156 values). Other parameters are tc = 10, r e ff =16 /un, T^A-o-SS/un = 0-2. 
(b) Reconstruction of W K using Eq. (11b), where slopes and intercepts are from the linear 
relationships in (a). All 211 wavelengths from 0.4 to 2.5 /xm are retrieved. 


be derived from these two variables and co Qx . Thus for 
estimating a broadband integral, a small number of ra- 
diative transfer calculations may be sufficient if high 
accuracy is not required. Note that W A = 1 - A a has 
been approximated as a function of all shortwave wave- 
lengths from 0.4 to 2.5 /urn using radiative transfer cal- 
culations at only three wavelengths (see Fig. 12b) with a 
bias of 0.006 out of 0.53 or 1.1%. However, if both the 
broadband /? A and T A are required, the biases will be 
bigger. 

c. Filling missing spectral data 

When high temporal and spectral resolution radiance 
measurements are provided by an aircraft spectrometer, 
some spectral data may be lost because of saturation 
or transmission problems (A. Vogelmann 2010, personal 
communication). In this case, if the SI relationship is 


established and the SI variables are found, the lost spec- 
tral data can be reconstructed using Eq. (4). 

d. Testing 3D radiative transfer codes 

The best way to test a new numerical radiative 
transfer model is to compare its results with a nontrivial 
exact solution of the radiative transfer equation. Un- 
fortunately, few such solutions are available, especially 
for 3D radiative transfer. To compare different atmo- 
spheric radiative transfer models, the Intercomparison 
of 3D Radiation Codes (I3RC; see http://i3rc.gsfc.nasa. 
gov and Cahalan et al. 2005) used some “consensus” 
results. To replace the unavailable “truth” for the Radi- 
ation Transfer Model Intercomparison (RAMI) in veg- 
etation, Pinty et al. (2001, 2004) introduced the “most 
credible solution.” However, both intercomparisons of 
radiation codes (I3RC and RAMI) mostly compared the 
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ox 


Fig. 11. RMS and bias [see Eq. (15) with n = 211] of W x for 
different approximations based on wqa thresholds. 

model results against each other. The SI theory suggests 
a more robust physics-based approach for testing radi- 
ative transfer codes. This approach states that the ratios 
of the calculated radiative quantities (radiances and/or 
fluxes) to the spectral single scattering albedo for three 
or more wavelengths are expected to be approximately 
linear in these radiative quantities. 

Note that the direct application of the SI relationships 
is limited only to cases with a purely absorbing surface. 
However, the relationships are an important part of the 
radiative transfer problem with general boundary con- 
ditions (see the appendix). More studies are needed to 
better understanding the limitations of applicability of 
the SI relationships to real atmospheres. 

7. Summary and discussion 

The solution of the radiative transfer equation de- 
pends on its boundary conditions (solar illumination 
and surface reflectance) and on three (independent) 
variables — optical depth r A , phase function P x , and sin- 
gle scattering albedo &> 0A . We assume that the functional 
dependence of the solution of the radiative transfer 
equation / A on {r A , P A } does not change with wave- 
length while the wavelength dependence of / A comes 
only from the oj oa spectra. In this case, we can formally 
write 

'a(V*>oa) = 7 a( t ’ P ; w oa)- ( 16 > 

Such a separation of variables is natural for radiative 
transfer in leaf canopies where scattering centers are 
much larger than the wavelength of solar radiation and 


dependence on t and P is determined by canopy struc- 
ture while dependence on a> 0A comes entirely from leaf 
physiology. 

For atmospheric radiative transfer this assumption is 
not met, in general (see Figs. 3a, b and 4a), since the size 
of scattering centers (air molecules, aerosol and cloud 
particles) is comparable to (or smaller than) the wave- 
length of solar radiation. 

However, in cloudy atmospheres assumption (16) can 
be met approximately for a large range of wavelengths. 
This paper shows that, at least for tq — 3 and wave- 
lengths that exclude strong water vapor absorbing 
spectral intervals, Eqs. (5a) and (5b) hold with variations 
less than 5 % (Table 1). Consequently, the linear spectral- 
invariant relationships are valid (see Figs. 6b and 7a, b); 
that is, the ratio of the radiative quantity (radiance or 
flux) to the total single scattering albedo is a linear func- 
tion of this quantity; its slope and intercept are wave- 
length independent. 

The slope of the linear function is an approximation 
to the recollision probability — the probability that a 
scattered photon will interact with the medium again. 
The intercept of the linear function quantifies the escape 
probability in a given direction. For hemispherically in- 
tegrated fluxes [Eqs. (8a) and (8b)], the sum of the slope 
and the intercept is equal to one [see Eqs. (7) and (12) and 
Fig. 6b]; this is not necessarily true for radiances (Figs. 
7a, b). 

The ratio of the absorbed radiation to the single 
scattering co albedo estimates the number of interac- 
tions; the reciprocal of the number of interactions is 
linearly related to the single scattering albedo [Eq. (9)]. 
As a function of wavelength, the total photon path 
length can be estimated [Eq. (10)] using explicit numer- 
ical radiative transfer models such as SBDART for ID 
and SHDOM for 3D calculations. 

The results presented in the paper can be directly 
applied to simple radiative transfer calculations if the 
same radiative quantity is required for different single 
scattering albedos (see Figs. 2a, b). In brief, we were able 
to single out one spectral variable — the single scattering 
albedo — that determines spectral behavior. If the single 
scattering albedo is known, it will help to fill spectral 
gaps in the spectral observations. In addition, the pos- 
sible applications of SI to climate and remote sensing 
research are summarized in section 6. 

If the spectral single scattering albedo is available, the 
slope of the linear relationship can be obtained directly 
from observations without any parameterizations and/or 
radiative transfer calculations. The product of the slope 
(the recollision probability) and the single scattering al- 
bedo approximates the maximum eigenvalue of the ra- 
diative transfer equation (see the appendix). As far as we 
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(b) wavelength, A (pm) 

FIG. 12. (a) Relationships of vs W A for A = 0.4, 1.6, and 2.1 pm. Other parameters are 

t c = 10, r e (f = 16, t a . a = 0^5,011 = 0.2. (b) Reconstruction of W A using Eq. (11a) or (7) with slope 
(asp) and intercept (as 1 - p) obtained from (a). Retrieved values for all 211 wavelengths from 
0.4 to 23 p.m are shown. 


are r.ware, this is the first method for estimating the 
maximum eigenvalue directly from atmospheric mea- 
surements. 

The concept of collision and escape probabilities was 
first introduced and developed in nuclear reactor physics 
(Bell and Glasstone 1970, 115-125). The vegetation 
community used this concept to interpret the SI re- 
lationships observed in radiative measurements and 
satellite data and to relate the wavelength-independent 
recollision and escape probabilities to canopy structure. 
In this paper the principles of SI have been applied to 
atmospheric radiative transfer. This provides a bridge 
between the methods developed in reactor physics, veg- 
etation canopies, and cloudy atmospheres since all three 
use the same radiative transfer equation. 

Finally, recently Marshak et al. (2009) reported a 
surprising SI relationship in shortwave spectrometer 
obser\ations taken by the Atmospheric Radiation 


Measurement program (ARM). The relationship sug- 
gests that the shortwave spectrum near cloud edges can 
be determined by a linear combination of zenith radi- 
ance spectra of the cloudy and clear regions. Chiu et al. 
(2010) confirmed these findings with intensive radiative 
transfer simulations of the different aerosols and clouds 
properties as well as the underlying surface types, and 
the finite field of view of the spectrometer. Although 
their calculations are performed for ID clouds, the first 
3D radiative transfer results suggest that the SI relation- 
ship discovered in shortwave spectrometer measure- 
ments is valid for the 3D simulation world. However, 
the slope and intercept depend on cloud structure and 
may be different from their ID counterparts. A clear 
physical (and mathematical) understanding of the ob- 
served and simulated SI behavior of zenith radiance 
around cloud edges is still missing. The current paper is 
the first step in this direction. 
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APPENDIX 

Eigenvalues of the Radiative Transfer Equation 
and Spectral Invariance 

a. Definition 

An eigenvalue of the 3D radiative transfer equation is 
a number y such that there exists a nontrivial function 
e(x, Q) that satisfies the equation 

y[fl • \e(x, 11) + o-(x)e(x, H)] 

o- a (iy - a)e(x,a') da', (ai> 

4n 

and zero boundary conditions. The 3D radiative transfer 
equation has a discrete set of eigenvalues y k and ei- 
genvectors e k {x, a), k = 0, 1, 2, Under some general 

conditions, the eigenvectors are mutually orthogonal. 
Solution of the radiative transfer equation can be ex- 
panded in eigenvectors. There is a positive eigenvalue 
y 0 that corresponds to a positive eigenvector e 0 (x. li). 
The eigenvalue yo is greater than the absolute magni- 
tudes of the remaining eigenvalues. It means that this 
eigenvalue and corresponding eigenvector dominates, all 
other eigenvalues and eigenvectors. The reciprocal of 
the positive eigenvalue y 0 describes the criticality con- 
dition in reactor physics — that is, a relation between 
l/y fl and the size of the assembly when more than one 
neutron is emitted per collision (Bell and Glasstone 
1970; Case and Zweifel 1967). Details of this theory can 
be found in Vladimirov (1963) and Knyazikhin et al. 
( 2011 ). 

The positive eigenvalue y 0 and corresponding eigen- 
vector e 0 have the following interpretation. Equation 
(Al) can be rewritten in terms of linear operators: dif- 
ferential operator L [the square brackets of Eq. (Al)] 
and integral operator S [the right-hand side of Eq. (Al)], 
as y 0 Le 0 = 5e 0 . Solving this equation for y 0 yields y 0 = 
L -1 Se o/e 0 . The operator L -1 5 inputs the source e 0 , sim- 
ulates the scattering event S and the photon free path 
L“\ and outputs the distribution of photons just before 
their next interaction. The ratio L~ l Se 0 /e 0 therefore can 


be treated as the probability that a scattered photon will 
recollide. 

b. Conditions for wavelength independency 

In cloudy atmospheres, the scattering coefficient <r SiA = 
i oo\aP , where ca 0A and P are the single scattering albedo 
and scattering phase function, respectively. The scat- 
tering operator S takes the form S - o) 0X S 0 , where So is 
determined by the extinction coefficient a and P. Equa- 
tion (Al) can be rewritten as pLe = S 0 e, where p = y/oj 0A 
is the probability that a scattered photon will recollide. 
If the extinction coefficient and scattering phase function 
do not depend on wavelength, L and 5 0 are wavelength 
independent. As a result, the recollision probability be- 
comes wavelength independent and does not depend on 
the incoming radiation either. 

c. Comments on Eq. (3) 

The 3D radiative transfer equation with a nonreflecting 
boundary can be transformed into the integral equation 

7 a = -oaVa + 2’ (A 2 ) 

where T 0 =L -1 S 0 ; the source term Q describes the 
distribution of uncollided radiation. The operator T 0 
has the same set of eigenvalues and eigenvectors as 
Eq. (Al), i.e., y *€* = a) 0x T 0 e k . If er does not depend on 
wavelength, Q becomes wavelength independent. In 
addition, if P does not vary with wavelength, T 0 becomes 
spectrally invariant. Solution of the radiative transfer 
equation can be expanded in successive order of scat- 
tering, which has a form of series similar to Eq. (3); that 
is, the terms in series (3) are powers of To that are pro- 
portional to poj™ K p m ~ x . More details on relationships 
among T 0 , recollision, and escape probabilities as well as 
their dependences on the scattering order m can be 
found in Huang et al. (2007). 

Note that the concept of spectral invariance is for- 
mulated for a medium bounded by a nonreflecting sur- 
face. The 3D radiative transfer problem with arbitrary 
boundary conditions can be expressed as a superposition 
of the solutions of some basic radiative transfer sub- 
problems with purely absorbing boundaries (Davis and 
Knyazikhin 2005) to which the spectral invariants are 
applicable. 

d. Number of interactions 

The mean number of photon interactions N x is an 
important variable characterizing the effect of multiple 
scattering on the atmosphere reflective and absorptive 
properties. This is just the integral of crl over spatial and 
angular variables; that is, 
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*A = IW 



c r(x)I A (a:, fl) fie d£l 


(A3) 


over a domain V. Multiplying Eq. (A2) by cr and inte- 
grating over spatial and angular variables yields 




(O 


0A 


IVaI 


iI 7 a 


-N, + L 


o- 


(A4) 


Here io = hv JV a ( X )Q(x,fl)dxdfl = 1 - T dir , where 
7dir is the directly transmitted radiation (section 4). 
Neglecting all terms in the expansion of the solution of 
the radiative transfer equation in eigenvectors except 
the dominant one e 0 , one gets j[ 7oe 0 i|/]|eol| *** P • Substituting 
it into Eq. (A4) and normalizing by i 0 , one obtains Eq. 
(9), where = NJi Q . Given the absorption can be 
estimated as (1 - w 0 a)N a . More details can be found in 
Panferov et al. (2001) and Knyazikhin et al. (2011). 


e. Energy conservation law 

Let A x and S x ~ R\ + T x — T di! be fractions of pho- 
tons from the incident beam that are absorbed and dif- 
fusely scattered (reflected or transmitted) out by the 
layer, respectively. The energy conservation law takes 
the form (Smolander and Stenberg 2005) 

+ *A + T k = A a + + r djr = 1. (A5) 

The diffusely scattered radiation is then 

= 0- ~ ^dir) ~ “ w Q>A*0 = J _ po} ^ ‘(V 

(A6) 

Equation (7) for the layer scattering W A = Sxlia directly 
follows from Eq. (A6). [Note that Ax in Eq. (A5) and 
A 0 x in Eq. (9) are different by a factor of t'o-] 
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